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We calculate how correlations in a Bose lattice gas grow during a finite speed ramp from the 
Mott to the Superfluid regime. We use an interacting doublon-holon model, applying a mean-field 
approach for implementing hard-core constraints between these degrees of freedom. Our solutions 
are valid in any dimension, and agree with experimental results and with DMRG calculations in 
one dimension. We find that the final energy density of the system drops quickly with increased 
ramp time for ramps shorter than one hopping time, JTramp fS 1- For longer ramps, the final energy 
density depends only weakly on ramp speed. We calculate the effects of inelastic light scattering 
during such ramps. 

PACS numbers: 


I. INTRODUCTION 


The dynamics of systems driven through a phase tran¬ 
sition are a source of rich physics [T]. The phenomenol¬ 
ogy is particularly interesting in zero-temperature sys¬ 
tems driven through a quantum phase transition [a [3]. 
In recent years, breakthrough experimental techniques 
in atomic physics have given us a direct probe of such 
transitions [US]. In this paper, we model a bosonic lat¬ 
tice system driven from a Mott insulator state to into 
the superfluid regime. We introduce a novel mean-field 
theory, building on commonly used doublon-holon mod¬ 
els |9]. We calculate how correlations develop during a 
lattice ramp through the phase transition. 

The phase diagram of bosonic lattice systems has been 
explored thoroughly [IDHIlj- In the strongly interacting 
regime, at commensurate filling, lattice bosons form an 
incompressible Mott insulator. Conversely, for weak in¬ 
teractions the ground state is a superfluid Bose-Einstein 
condensate with long range order. When the system be¬ 
gins in a Mott insulator state and interactions are turned 
off, correlations grow as quasiparticles propagate across 
the system nai^. 

The Mott and superfluid phases can be approximated 
by distinct mean-field quasiparticle models. The excita¬ 
tions in the superfluid phase are well described by Bogoli- 
ubov quasiparticles made up of superpositions of parti¬ 
cles and holes m- In the Mott insulator regime, on-site 
number fluctuations are small and the occupation of each 
site can be truncated to a small number of possibilities 
[HI, the “doublon-holon” model. At strong coupling, the 
doublons and holons can be approximated as noninter¬ 
acting bosons. These two descriptions are incompatible, 
making it challenge to model the dynamics across the 
phase boundary. 

Previous work has produced partial understanding of 
this transition[T8l[T9]. Product state methods such as the 
Gutzwiller ansatz cannot calculate correlations [201 El] . 
Other approaches have included calculations on small 


lattices [55], field theory calculations for large particle 
density pHVEi)] and various numerical techniques, which 
work well in one dimension but are otherwise more lim¬ 
ited |57H21]. There has also been significant work on 
sudden quenches EHIISI] Here, we provide an analytical 
model that is particularly suitable for the small mean oc¬ 
cupation numbers common in atomic experiments, pro¬ 
vides access to coherence data, and is applicable in any 
number of dimensions. 


II. MODEL 


We perform our calculation within an approxi¬ 
mate doublon-holon model. We restrict the state 
of each site i to the subspace of occupations |i) G 
{|h -|- 1), |h), \n — I)}, where n is the median number of 
particles per site. The system can then be thought of in 
terms of a mean-occupation background and hard-core 
quasiparticle excitations of “holons” (an h — I occupa¬ 
tion) and “doublons” (h -I- 1 occupation). The annihi¬ 
lation operators at site i for these quasiparticles are de¬ 
fined by di\n)^ = di\n-l)- = hi\n + l)- = hi\n)- = 0, 
di\n + I)j = hi\n- I)^ = \n)^. 

Under this approximation, the Hamiltonian is 


k 




dpdf^ -f hih}^ 


ti 


£k {d\dk - hlhk^ + JnSk (dkh-k + ■ 


( 1 ) 


Here, dk = e-^^'^'di, summing over all sites i, and 

similar for hk, while Sk = —2 cos(fe • A), summing 

over lattice basis vectors, A = Ax, Ay, Az in three di¬ 
mensions, or a subset of those in lower dimensions. These 
represent a cubic lattice with lattice constant A. U and 
J are the in teraction and hopping strength, respectively, 
and h = n{n + 1). Ns is the number of sites in the 
lattice. 
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The doublon-holon model is an approximation for the 
single-band Bose-Hubbard model [TOl [32] . It is most ac¬ 
curate in the low-temperature, strongly-interacting limit, 
as the energy of a state increases quadratically with the 
deviation from the mean particle number. However, for 
low occupation numbers h, it can be a good approxima¬ 
tion in the weakly-interacting limit as well. In a non¬ 
interacting superfluid gas with h = 1, the probability of 
finding more than two particles on a given site is less than 
10%. We do all our calculations in this regime, taking, 

{hi) =h + - (h\hi^ =h = l. 

We calculate the time evolution of the two-point corre¬ 
lation functions, (dldk) ~ Idkh-k) using 


the Heisenberg equation, 4(X) = i H,X 


The hard¬ 


core constraints for di,hi imply nontrivial commutation 
relations, and the resulting equations of motion involve 
four-point correlation functions such as 


— (^^ph_p_qh—k—qd}^. ( 2 ) 

We can characterize Ck,p,q by writing it in the form 



This 

equation defines 

the function 

^p^k,q- 

We 

make a mean-field 

approximation. 

tak- 

ing 

^p,k,q ~ -^(h^-p-qh- 

p—q''j (^d^pdp'j, 

where 

n‘^ = 

l^Efc(44) is the 

doublon density. 

This 


approximation enforces the hard-core constraint 
J2k ^k,p,q — 0, and becomes exact in the deep Mott 
regime. We make similar approximation for the other 
four-point correlation functions, as described in detail in 
Appendix 

We arrive at a closed set of non-linear, coupled dif¬ 
ferential equations that we numerically integrate to find 
all quasiparticle two-point correlation functions at any 
time. From these we can easily extract the correlation 

functions for real particles, (dlaj) and la\au\ 



FIG. 1: (Color Online) Equilibrium properties of the hard¬ 
core doublon-holon model discussed in the text for a cubic 
lattice with mean filling h = 1. Shown as a function of the 
interaction strength U/J, above: the equilibrium condensate 
fraction, below: the average energy per particle in the ground 
state. 



FIG. 2: (Color Online) Evolution of the momentum density 
distribution function as the interaction strength is 

slowly ramped down, U = Ui{Uf in a one-dimensional 

lattice. Here Ui = 47 J, H/ = 2J, Jtt = 2. 


III. EQUILIBRIUM STATE 


IV. INTERACTION RAMPS 


We find the equilibrium state under this model by min¬ 
imizing the expectation value of the Hamiltonian of 

Eq. (1) while requiring f^(d\dk^ = f^(dkh-k^ = 0. As 
seen m Fig. we hnd a phase transition at a critical 
value of Uc/J = 10.4,21.8,33.4 in one-, two- and three- 
dimensions. These are similar to the standard mean- 
field values of Uc/J = 11.6,23.2,34.8 [TT1|33| and some¬ 
what higher than numerically calculated values Ud J = 
3.6,16.9,29.3 


We use the model above to explore the behavior of a 
gas subject to a non-adiabatic ramp of the interaction 
through the phase transition. We perform an interaction 
ramp of the form 


V = Ui{UdU/fl^U (4) 

where the ground state of the system is a Mott insulator 
for U = Ui and superfluid for U = Uf. The time scale 
Tr sets the speed of the ramp. This form approximates 
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the relation U/ J in an optical lattice experiment if the 
scattering length is fixed and the lattice depth is ramped 
down [T^ . 

We initialize the system in the ground state at the 
initial lattice depth, in the Mott regime, and perform a 
finite-element time integration of the evolution equations 
as the interaction strength is reduced. We calculate the 
momentum space density throughout this evolution for 
various values of Xr. Figure shows the behavior for a 
typical ramp, with Jtt = 2. We have full access to all 
two-point observables at any time along the ramp. 

We first characterize the behavior of the system at 
the end of the ramp. We define an effective correlation 
length, by comparing correlations in the system to the 
form We calculate ^ by fitting to 

the width of the momentum distribution, as defined by 
the first moment, yielding 


i 

A 


-1/log 



Though it is infinite for an equilibrium superfluid system, 
^ remains finite at any finite time for a system that is not 
initially superfluid m- 

Figure shows the effective correlation length at the 
end of the ramp for varying ramp times. In one dimen¬ 
sion our calculation agrees well with the result of an ex¬ 
act diagonalization of a small lattice. The discrepancy is 
consistent with the finite-size effects in the exact diago¬ 
nalization. Our results also agree with the experimental 
results of [40] . For slow ramps we see that the correlation 
length is somewhat smaller in lower dimensions. 



V. FINAL ENERGY DENSITY 


After the ramp has ended, the system continues to 
evolve, and the correlation length continues to grow. 
However, the energy of the system is now conserved. At 
long times after the ramp we expect the state of the sys¬ 
tem to resemble a thermal state at a temperature de¬ 
termined by the energy density U = , 

where \H) is the energy of the new ground state of 
the system. 

We plot U as a, function of the ramp time Tr in Fig. [^ 
For ramp times much shorter than the hopping time 
scale, jTr < 0.2, the final energy density varies slowly 
with Tf. Such ramps are indistinguishable from instan¬ 
taneous quenches, and the final state of the system, if 
allowed to equilibrate, would be similar for any Tr in this 
regime. For Jr^. > 0.2, the system’s energy depends more 
strongly on the length of the ramp. 

Figure]^ also shows the critical energy density Uc cor¬ 
responding to the energy density of of a Bose lattice gas 
with U = 17/ at the critical temperature of the superfluid- 
normal gas phase transition [371138]. We expect a gas 
with U > Uc to equilibrate to a normal-gas state with 
finite correlation length while a gas at U <Uc would 
equilibrate to a superfluid state with long-range order. 
In two dimensions, we expect short ramps, Jr^ 0.6, 
to lead to a normal state, while longer ramps lead to a 
superfluid gas. In three dimensions, the energy density 
is always below lAc, even for an instantaneous quench. In 
one dimension there is no condensed phase. 



FIG. 3: (Color Online) Effective correlation length ^ (see 
Eq. ©), normalized by the lattice constant A, at the 
end of a ramp of the interaction strength of the form 
U = C/i(17//17„)‘''"'-. Here Ui = 47 J, Uf = 2J. The red dots 
are the result of an exact diagonalization calculation for an 
11-site one-dimensional lattice. 


FIG. 4: (Color Online) The energy density U following an 
interaction ramp of the form U = Ui{Uf . Here 

Ui = 47J, Uf = 2J. Horizontal lines show the energy den¬ 
sity, UajJ = 0,2.1, 5.1, at the superfluid critical temperature 
Tc/J = 0,1.7, 5.9, for U = Uf = 2J, in one, two and three 
dimensions [371HH]. 
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VI. DECOHERENCE 


In an ideal, closed, quantum system, all evolu¬ 
tion is unitary. The final energy of the system rises 
monotonously with the rate of the ramp in such systems. 
Conversely, any real system suffers from heating, atom 
loss and other impacts from the environment. As a re¬ 
sult, experimental dynamic systems always face a com¬ 
petition between the system’s reaction time and external 
processes. 

The physics of such decoherence has been explored 
in detail [41H46| . Here, we return to a mechanism we 
have previously used to described the effect of density 
measurement by light scattering m- The same formal¬ 
ism describes inelastic light scattering, where an external 
photon scatters off of a trapped atom. This is one of the 
major sources of decoherence in atomic experiments. 

As in SZ], we neglect out of band effects, which cause 
particle loss. We focus on in-band scattering, which 
would directly decrease the coherence of the remaining 
gas and reduce the correlation length measured above. 
In an ensemble description, this leads to a nonunitary 
evolution term of the form 






— {dkh-k) = -i 


dk^ — k 5 H 


{(S'j^dk J - rid 
- lidkhk 


( 6 ) 


where 7 is proportional to the frequency of light scatter¬ 
ing per site. 

We calculate the effect of this decoherence on the be¬ 
havior of the correlation length as shown in Fig. 

As expected, no effect is seen at time scales shorter than 
1 / 7 , but at longer time scales, inelastic processes cause 
the correlation length to decay. The overall effect is sim¬ 
ilar to experimental observations in [4(1] . 



FIG. 5: (Color Online) Effective correlation length at the 
end of a ramp of the interaction strength, in a system cou¬ 
pled to the environment in the form shown in Eq. §. Here 
Uf = 47J, t/i = 2J, in a three-dimensional cubic lattice. In 
an optical lattice setup, the rate of inelastic light scattering 
events changes with lattice depth similarly to the interaction 

strength [sms]- 
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VII. OUTLOOK 

The physics of ultracold atomic systems involves multi¬ 
ple energy scales. In driven experimental systems, these 
include the relaxation time of the system, the driving 
time scale and the rate of decoherence imposed by inter¬ 
action with the environment. Here, we have quantified 
the effect of the quench rate in Bose-Hubbard systems 
crossing the phase boundary. We find that there are two 
regimes. For sweeps which are much shorter than the typ¬ 
ical hopping time, Jr^ < 0 . 2 , the ramp time has no effect 
on the final state and the ramp is indistinguishable from 
an instantaneous quench. For longer ramps, the final en¬ 
ergy density of the state and therefore its correlations at 
equilibrium, depend on the length of the ramp. In two 
dimensions, shorter ramps lead to a normal gas state, 
while longer ramps lead to a superfluid state. We have 
also demonstrated that inelastic light scattering can be 
quite destructive on longer time scale, underscoring the 
usefulness of shorter experimental runs. 
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Appendix A: Detailed Derivation of the Approximations in the Hard-Core Doublon-Holon Model 


1. Underlying Model 


We perform our calculation within an approximate “doublon-holon” model. The state of each site i can be given 
in terms of a spinor in the allowed occupation states \n + 1)^, \n — 1)^ where n is the median number of particles 
per site. We define the quasiparticle annihilation operators, di = |n)(n + l|j, hi = |n)(n — l|j. 

Under this approximation, the Hamiltonian is 


H = - 


J ^ {{n + l)d\dj + nh\hj + y/n(n + l){dihj + + h. c.^ + ^ + ^1^*) 


(*j> 


f f + J\l\ {d\dk + h^Jif^j + ^ek(dldk - h\hj^ + Jhek[dkh-k + 


(Al) 


Here, 




1 pjfe ri V 

. nr / r j 




Ah-Vi 


(A2) 


summing over all sites *, and 


£k = —2 cos(fc • A) 

A 


(A3) 


summing over lattice vectors. We perform our calculation on a cubic lattice with lattice spacing A, A = Ax, Ay, A£ 
in three dimensions or a subset of those in lower dimensions. U and J are the interaction and hopping strength, 
respectively, and h = \/n{h + 1). Ng is the number of sites in the lattice. 

We do all our calculations for a density of one particle per site, {hi) =h + (S'idi^ — =n = l. 

The hard-core constraints on the operators d, h translate into non-trivial commutation relations. 


dk,dl 


= 5k,q - 2nf_t. - 


^q—k 


q—k^ 


r}'^ h 

fiq 


hkAl — 5k,q h„_k 2h„_k, 


— ^q—k 
= ^k-q 


^k 5 


hk, hq 


dk') hq 


= 0 , 


(A4) 


where we define the quasiparticle density operators 




''d]di 




— ik-‘ 


‘h]h. 


_ 


= iiE' 


— ik‘ 


’''hU,. 


(A5) 


We write hd,h = the density of doublons and holons, respectively. In the Mott equilibrium limit, the operators 
in Eq. (A5) can be neglected and the quasiparticles can be treated as noninteracting bosons. This is not true in the 
superfluid regime. 


2. Equations of Motion 


Equations of motion can be derived from the Hamiltonian, Eq. (Al), via the Heisenberg equation. 


_d 

dt 



= i 


X,H 


(A6) 
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We focus on the two-point observables, 


_d 

dt 



dt 


(hlkh-k) 


d 

dt 


dkh-h 


iJnek{(dkh_k^ - (d\h'_j^^ 




(\/i++ m) 

{dl{^<-q 

+ n 

t,)4) 

— iJn 


+ (\/i + im ■ 

- 2 k) (hlql'-k- 

-gdfc) -h.c. 


Q 

_+((2n^. 

+ f^k-q 

)d-qdkj + 


q-kdqdk^ 

iJfiSk 

(^(^dkh 

-.) - (4'.h)) 






(\/i + ^~ m) 

(4,(2*;- 

k + 

n'^q-k)d-k) 

— iJh 


+ (\/l + 

+ 2k)(4^V 



Q 

_+((2<-fe 

+ K-k 

j dqh-k^ + 

(*; 

q_kh-qh-k^ 

-iu(^ 

dkh-k^ 

) - iJn£k{l 

- ind - 

?>nh) 



— iJh 


dk,q ~ 1^) 


ihi^dqh-q^ 

)+ 

{dldq'j + (hlqh- 


(\/l + 4^ + m) + (dqV-q-kdkYj 

+ (\/l + 4^-^l + 2hJ-fc)^-fc) + {h-ql'lg_kh-k)) 

+ {hlg(2n‘^_^ + n'^_M_k) + 

+ (4(«Lg + ‘2nt-q)dk) + (h^_qV-q-kdk) 

- ( 2 h^_fc -b -b 2 h^_ J - i -b 


-b * Jh f 


Here h. c. stands for the Hermitian conjugate. 


(A7) 


3. Hard-Core Coherent Approximation 

To perform the time evolution, we must make approximations for the quartic terms, such as 
Ck = ^ ] ^q\ d\nf._^dk) = 'y ^ y ] £qCk,q,p- 


p,q 


p,g 


These terms can be written out as 

^k ~ 7w y , ^k{djJl_k_jT--k-pdk) + 7w ^ ^, ^q{d\h_Jl-kdk) 

p 9 

~Wt,^ki^djJp_jJi—kd}^-\-j^ y ^ Eq(^dl^h_p_^h—k—qdk'^- 


p,g 
q^k,p^0 


(A8) 


(A9) 


The first two sums on the right hand side add up coherently, and we expect them to dominate. The third term is 
inversely proportional to the system size, and is therefore negligible. For bosonic operators, one may expect the final 
sum to add up incoherently, as in the Hartree-Fock-Bogoliubov approximation m, suggesting the form 


ci^ci = {^Y.ChU_;h.k- 


\ek{didk)+ 


h-kdk 


(AlO) 
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This intuition fails in the hard-core case. This can be seen by summing over the momenta, 

= £q(dlhlp_gh_p_kdkj = e"^'"'Sq{dlhlp_ghidi^ = 0 


p,q,k 




k \ P / \ k ) 

To account for the hard core constraints, we formally write 

cip^q = idXhip_;h_p_A) = 4.,(/iV9^-p-fc)44+ dp^^id\hi^l^h.uk) - did 




where this equation defines ak,p,q- We approximate this function with the hard core constraint in mind, 

1 


k^k^p^q 






h — T) — 


-p-q'^-p-q 


SO that J2k ^kpq — We then find 


4 = 


dlk)^ 

){kdk^ +£o'n*{dkh^k 


where 


nd = {nk ^d=i;Et(44) d = 


h_f,h 


qll_p_gh-p-f^dk 


^q,k(^h^—p—k^~P~k^(^^kh—k^ H“ ^p ,0 fc ^ (^h^—p—qh—p—q^(^dqh—q^(^^f^dk 


h—qd^_„dp—kdk 


q^p-q^P 

q,k{dl_kUp-k 


5n.k( dkuk-k')( dkh-k) + dj,-a.k( dad-a') ( kdk) k ud {dqd-q) (dp-qdp-q 


kkp-\-qdp—kd}^ ~ Sp^ol^kk^^kd—k'^ k+q^kl^qd—q'^(^k^dk 


- kk( dnh^ 


Ns rid\ ^ ^ 


d\.k 


''p+q^P+q 


Applying these approximations to Eq. (A7), we find 


d 

dt 


dldk'j = iJhek(^(^dkh_k'j - (^^kd-k")^ ~ 


^(^Sk^d(^dkh_f^^ SqT] (dik)) 

+2^0 + kfid* {dkd-kj + ^d{dkh-kj^ 


dp,dk 


— h. c. 


dt 


kh-k^ = -iu(^kh-k^ - iJnek(l - Grid + 9(n^ - -I- 6|77|^^ 


— 2iJh 


2iJh 


\J^ + 4^2 {J\dq^-‘^ vj + (^kk'j ^d 


1 + XST 


{i(^£knd(dkh-k^ - So'nl^dlk^') + 2£o4<^4^-/c^) 
-l-3(efcnd - eo^d){dlk) + 3£oT]*(kh-k 


with (h^_f.h-k'^ = (^dlky 


(All) 


(A12) 

(A13) 

(AM) 

(A15) 


We make similar approximations for the other terms, 

(^kqkp-qk-kk^ ~ Sq^k(dp-kk-k'^(^kk'^ + dp-q^klyik'^ikk'^ — wi:rdi(^kk'j(dp-qk~'y(d^kk^ 
Jkqdp_^ljhp—kdk'^ ~ dpfi(jkqky(^dkh—k''^ dpJ^q^kik—qh—(y(^kdk'j ~ ^(jkqdq'^ (ylt+qk+q'^ (y^kk'^ 


(A16) 


(A17) 
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